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Numerical algorithms based on variational and symplectic integrators exhibit special features that 
make them promising candidates for application to general relativity and other constrained Hamil- 
tonian systems. This paper lays part of the foundation for such applications. The midpoint rule 
for Hamilton's equations is examined from the perspectives of variational and symplectic integra- 
tors. It is shown that the midpoint rule preserves the symplectic form, conserves Noether charges, 
and exhibits excellent long-term energy behavior. The energy behavior is explained by the result, 
shown here, that the midpoint rule exactly conserves a phase space function that is close to the 
Hamiltonian. The presentation includes several examples. 
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I. INTRODUCTION 

This is the first in a series of papers that explore the 
possible advantages of using variational and symplectic 
numerical integration techniques for constrained Hamil- 
tonian systems. A constrained Hamiltonian system is the 
Hamiltonian formulation of a gauge theory yj. For such 
a theory the canonical momenta, defined by the deriva- 
tives of the Lagrangian with respect to the velocities, are 
not invertible for the velocities as functions of the coor- 
dinates and momenta. As Dirac showed |2(, this implies 
the presence of constraints among the coordinates and 
momenta. The constraints are the canonical generators 
of the gauge symmetry. They appear in the action as 
part of the Hamiltonian, accompanied by undetermined 
multipliers. 

Constrained Hamiltonian systems are common in 
physics. Examples include electrodynamics, Yang-Mills 
theories, string theory, and general relativity. The nu- 
merical integration of Maxwell's equations for electrody- 
namics has been well studied. For example, with the 
finite difference time domain (FDTD) method, the elec- 
tric and magnetic fields are evolved using discrete forms 
of Ampere's and Faraday's laws Q. The FDTD dis- 
cretization automatically preserves the two Gauss's law 
constraints in the source free case. Yang-Mills and string 
theories are primarily used to describe elementary quan- 
tum systems, so for these theories the classical solutions 
do not play a critical role. Correspondingly, numerical 
methods for evolving the classical Yang-Mills fields and 
classical strings have not been thoroughly explored. 

The most challenging example of a constrained Hamil- 
tonian system, and the one that serves as my primary mo- 
tivation for this investigation, is general relativity. There 
is currently a great deal of interest in developing nu- 
merical methods for solving Einstein's equations. This 
interest is driven by recent advances on the experimen- 
tal front. A number of ground-based gravitational wave 
detectors are in operation today, and during the next 
decade some of these instruments will reach the level of 
sensitivity needed to detect black hole collisions. The 
LISA project is a joint effort between NASA and ESA, 
with the goal of placing a gravitational wave detector in 



solar orbit. The LISA detector will be capable of sens- 
ing, among other sources, collisions between the super- 
massive black holes that reside at the centers of galaxies. 
To maximize the scientific payoff of these instruments 
we need a theoretical understanding of the gravitational- 
wave signals produced by black hole collisions and other 
astrophysical phenomena. The only known method for 
predicting the gravitational wave signature of colliding 
black holes is through numerical simulation. 

Numerical relativity is not a mature field. Researchers 
have spend much time and effort in developing numerical 
relativity codes, but the complexity of the Einstein equa- 
tions coupled with the topological issues that arise when 
modeling black holes have made progress slow. Current 
codes can succeed in simulating at most about one orbit 
of a binary black hole system before errors completely 
spoil the results 4-]. The main difficulty appears to be 
the presence of "constraint violating modes" [1, IE S • 
These are solutions of the Einstein evolution equations 
that are unphysical in that they do not respect the con- 
straints. Although the evolution equations preserve the 
constraints at an analytical level, numerical errors in- 
evitably excite these constraint violating modes. Some 
of these modes grow exponentially fast and spoil the nu- 
merical results. What is needed for numerical relativity 
is an algorithm that will keep the constraints satisfied, 
or nearly satisfied, during the course of the evolution. It 
might be possible to develop a scheme like the FDTD 
method of electrodynamics, but the complexity and non- 
linearity of the Einstein equations makes this a difficult 
task. Some progress along these lines has been made by 
Meier 

In this paper I begin to explore a different route for 
keeping the constraints satisfied for general relativity 
and other constrained Hamiltonian systems. The idea 
is based on the use of variational integrators (VI). In 
the traditional approach to numerical modeling by finite 
differences, the continuum equations of motion are dis- 
cretized by replacing derivatives with finite difference ap- 
proximations. In the VI approach we first discretize the 
action, then derive the discrete equations of motion by 
cxtrcmizing the action. This approach was pioneered by 
a number of researchers beginning in the 1960's; for a 
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brief historical overview, see Ref . |l(J . Variational inte- 
grators have been developed further in recent years by 
Marsden and collaborators [ToL fTTj| . 

One of the key properties of variational integrators 
is that they are symplectic. This means that the dis- 
crete time evolution defined by the VI equations auto- 
matically conserves a symplectic form. The subject of 
symplectic integrators is well-developed; for an overview, 
see Ref. jl^. Variational integrators also conserve the 
charges associated with symmetries via Noether's the- 
orem. For our present purposes, the most interesting 
characteristic of variational and symplectic integrators is 
their behavior regarding energy. Although these integra- 
tors do not typically conserve energy, they exhibit excel- 
lent long-time energy behavior. For other integrators the 
energy errors typically increase unboundedly in time. For 
variational and symplectic integrators the energy error is 
typically bounded in time. 

There are various ways that one can develop a vari- 
ational integrator for constrained Hamiltonian systems. 
For example, one can extremize the action while keeping 
the undetermined multipliers fixed. In that case the con- 
straints will not remain zero under the discrete time evo- 
lution. But there is reason to believe that in many cases 
the constraint errors, like energy, will remain bounded in 
time [T5 | . Another option is to extremize the action with 
respect to the undetermined multipliers as well as the 
canonical coordinates and momenta. This is the most 
attractive approach from a number of perspectives. In 
this case the discrete constraints are imposed as equa- 
tions of motion at each timestep, so they are guaranteed 
to hold under the discrete time evolution. The trade off 
is that the undetermined multipliers of the continuum 
theory are actually determined by the discrete equations 
of motion. 

In general relativity the constraints cannot be solved 
for the multipliers unless the coordinates and momenta 
are chosen appropriately. The traditional choice of 
canonical coordinates [1-4| . the spatial metric, leads to 
generically ill-defined equations for the multipliers. Re- 
cently Pfeiffer and York 0, ^| have rewritten the con- 
straints using the conformal metric and the trace of 
the extrinsic curvature as coordinates. They show that 
the resulting equations for the multipliers are generically 
well-defined. In Ref. I rewrote the action and evo- 
lution equations in terms of these new coordinates and 
their conjugate momenta. This is one form of the ac- 
tion that is suitable for the development of a variational 
integrator for general relativity. 

The essential idea of using a discrete action to define a 
set of discrete equations of motion that both respect the 
constraints and determine the multipliers has also been 
studied in the context of gen eral relativity by Di Bartolo, 
Gambini and Pullin [H [H HJ El El ■ They refer 
to their approach as "consistent discretization" . They 
discuss consistent discretization in the context of numer- 
ical relativity, and also as a route toward quantization. 
There are a number of technical differences between the 



works of Di Bartolo, Gambini and Pullin and the results 
presented in this and the following papers. The most im- 
portant difference between my approach and theirs is a 
difference in techniques used to generate the equations 
of motion. I extremize the discrete action directly while 
Di Bartolo et al. identify the discrete Lagrangian as the 
generator of a Type 1 canonical transformation. With 
direct extremization we obtain useful information about 
the system encoded in the endpoints of the varied action. 
This is the key to proving the important properties of the 
variational integrator including symplecticity, Noether's 
theorem, and the good long-time behavior of energy. 

In this first paper I focus on simple mechanical sys- 
tems with no constraints. This is a rich subject that has 
been explored rather thoroughly, in mathematically pre- 
cise language, by Marsden et al. 0, 0] . The purpose 
of this paper is to present the key results on variational 
integrators in the context of a particular discretization 
using language familiar to most physicists. The partic- 
ular discretization of the action considered here leads to 
the midpoint rule applied to Hamilton's equations. The 
midpoint rule is an old, familiar numerical algorithm. It 
is presented here in a new, perhaps unfamiliar light as a 
variational-symplectic integrator. This new perspective 
allows us to derive and to understand the characteristic 
features of this integrator on a rather deep level. 

In the next section I review the derivation of Hamil- 
ton's equations from the action expressed in Hamiltonian 
form. In Sec. IIIII I discretize the action and derive the 
VI equations from its extremum. In Sec. II VI I show that 
the variational integrator is symplectic, and Noether's 
theorem holds. I also show that the VI equations can 
be written as the midpoint rule applied to Hamilton's 
equations. Section contains a discussion of energy. 
There, it is shown that the energy is well behaved be- 
cause the VI equations exactly conserve the value of a 
phase space function that is close, in a sense to be dis- 
cussed, to the Hamiltonian. Several examples are given 
in Sec. IVII These examples explore the energy behavior 
and the convergence properties of the midpoint rule as a 
variational integrator. 

My goal is to investigate variational and symplectic 
integration techniques for constrained Hamiltonian sys- 
tems. In the next paper in this series [13|, I will apply 
these techniques to a class of simple constrained Hamil- 
tonian systems, namely, parametrized Hamiltonian me- 
chanics. These are ordinary Hamiltonian systems with 
the coordinates, momenta, and time expressed as func- 
tions of an arbitrary parameter. The theory is invariant 
under changes of the parameter, and this gauge invari- 
ance gives rise to a constraint that enforces conservation 
of energy. In future papers I will apply VI techniques to 
field theories with gauge symmetries. In canonical form 
these theories are described as constrained Hamiltonian 
systems with constraints that are local functions in space. 
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II. CONTINUUM MECHANICS 

Let the index a label pairs of canonically conjugate 
dynamical variables x a and p a . The action is a functional 
of x a (t) and p a (t), given by 



S[p, x] 



dt \p a x a - H(p, X, t)] 



(1) 



Here, H (p, x, t) is the Hamiltonian and t is physical time. 
The dot denotes differentiation with respect t. The sum- 
mation convention is used for repeated indices, so the 
expression p a x a includes an implied sum over a. 
Variation of the action yields 



SS\p,x) 



t" 



dt 



PaSXa 



- tt-) 5x a 

t" 



(2) 



With the coordinates x a fixed at the initial and fi- 
nal times, t' and t" , the endpoint terms in SS vanish. 
Then the condition that the action should be stationary, 
6S = 0, yields 



Pa 



dH 

dp a 

dH 

8x a 



(3a) 
(3b) 



These are the familiar Hamilton's equations. An imme- 
diate consequence of these equations is that the Hamil- 
tonian function H(p,x 1 t) satisfies 



H 



dH 
~dt 



(4) 



If H has no explicit t dependence, then H = 0. In this 
case H, the energy, is a constant of the motion. 



III. DISCRETE MECHANICS 

Let us divide the time interval between t' and t" into 
N equal subintervals, or "zones", labeled n = 1, . . . , N. 
These zones are separated by nodes, which are labeled 
n = 0, . . . , N. As seen in Fig. the endpoints of zone 
n are nodes n — 1 and n. The expression t n denotes the 
time at the nth node. Likewise, x™ and p™ denote the 
coordinates and momenta at the nth node. The timestep 
is At = t n — t" _1 . In this paper I consider the following 
second order accurate discretization of the action QJ: 



JV 



S\p,x] = J2 At 



Axl 
At 



(5) 







\ 



time 



zone 1 



zone 2 



FIG. 1: Discretization in time. The nodes are labeled n = 
0, . . . , N and the zones (or time intervals) are labeled n = 
1, . . . ,N. The coordinates and momenta are node centered, 
the Hamiltonian function is zone centered. 



The A notation and the underlined index notation are 
employed repeatedly below; they are defined by 



Ax^ = 



-xl 



(6a) 
(6b) 



These operations commute; that is, 



Aj 



L )/2- 



It will also prove useful to denote the value of the 
Hamiltonian in the nth zone by 



H n 



H^x^) 



(7) 



That is, we view t, x a , and p a as node centered in time 
and H as zone centered in time. [See Fig. {Q.] Then 
equation expresses the fact that, to second order ac- 
curacy, the zone centered values of t, x a , and p a that 
appear in H n can be obtained from the averages of the 
neighboring node centered values. 

The discrete "Lagrangian" , that is, the term in square 
brackets in Eq. has truncation errors that scale like 
0(At 2 ). The discrete action is a sum over N ~ 1/At 
terms, each having errors of order 0(Ai 3 ). It follows 
that the error in S typically scales like 0(At 2 ). Thus the 
action JSJ is second order accurate. Note that Eq. (J5J 
is not the only possible second order discretization of 
the action. For pedagogical purposes, I have chosen to 
restrict considerations in this paper to the discrete action 
(j3J). Other discretizations, including some with higher 
order accuracy, will be discussed elsewhere [l3|. 



Note that the momentum variables appear in the ac- 
tion JSJ only in the combination p^t = (p™ +p" _1 )/2. 
This combination represents the zone centered momen- 
tum, accurate to second order. Let us set this observation 
aside for the moment and treat the action as a function 
of all node-centered coordinates and momenta, x™ and 
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p™ for n — 0, . . . , N. The variation of S is 
w-i 



SS = J2 At 



77=1 

N-l 

77=1 
1 

+ 2 



— 
At 



a^ 



71+1 



At 



\dXaJ 



8x n „ 



+- 



1 



Ax 



3p a 



JV 



At 



N 



' , i /saV . 



fa! 



„ 1 fdH\" , 



<5.r 



A' 



(8) 



Here and below we treat the derivatives of H(p, x, t), like 
H itself, as zone centered quantities. Recall the notation 
defined in Eqs. © and (JJJ. For any zone-centered func- 
tion F(p, x, t) of the canonical variables and time, we 
have F-±i ee [F n+1 + F n }/2 = [F(pS±l, x-±I, i-±i) + 
F{p rL , x— , i— )]/2. These notational rules apply to the 
derivatives of the Hamiltonian that appear in SS. 

If we fix the coordinates at the endpoints, x° a and x^, 
then the condition that the discrete action should be ex- 
tremized is 



AxP^~ 

At 

At 
Ax 1 



71 + 1 



71+1 



, n = 1, . . . , N-l , 



At 

Ax N 
At 



^] • »= I V I 



( — ) 

\dx a J 

dpa] 

dH_^ N 

dp a 



(9a) 
(9b) 
(9c) 
(9d) 



These equations are redundant. For example, Eq. (|9d|) 
can be derived from Eqs. (|5Jl,c). This redundancy is a 
result of the fact that the action does not depend on 
the node-centered momenta p™ independently, but only 
on the zone-centered combinations pit. We can combine 
equations 101, c,d) into a single expression and write the 
equations of motion @ as 



77 + 1 



Aa;™ +1 fdH 
At ~ \dpZ 

\dx a ) 



A n+1 
Apg 

At 



, n = 0,...,N-1 , (10a) 
n = I,..., N-l , (10b) 



The equations of motion in this form can be obtained di- 
rectly from the action (JjjJ by extremizing with respect to 
the node-centered coordinates x" and the zone-centered 



momenta pt. They are a discrete form of Hamilton's 
equations ©. 

The equations of motion l|10|) constitute the varia- 
tional integrator defined by the discrete action (0 ■ Since 
they are derived from a variational principle, these equa- 
tions naturally define a boundary value problem in which 
the freely chosen data are divided between the end- 
points in time. Thus, given the boundary data x° a and 



Eqs. ((Tn|) determine the coordinates x™ for 



1, . . . , N — 1 and momenta pit for n — 1, . . . , N. We 
can add boundary terms to the action to change the per- 
mitted boundary conditions. However, in practice, our 
primary interest is not in any of these boundary value 
problems. Rather, we are interested in solving an initial 
value problem. Thus, we are faced with the task of rein- 
terpreting the equations of motion in such a way that 
initial data can be posed and evolved into the future. 

It is not difficult to reinterpret the variational integra- 
tor (|1C)|> as an initial value problem. One possibility is 
to choose values for the coordinates at the initial time t° 
and values for the momenta at the half timestep t—\ that 
is, we choose x® and p\. Then Eq. (jlOaf) with n = can 
be solved for x\. This completes the determination of 
data at "levels" n = and 1. Alternatively, we can gen- 
erate data at levels n — and 1 by specifying x° a and x\, 
then solving Eq. (|10afl with n = for p^. Once the data 
at levels and 1 have been found, we can solve Eqs. i(TT)|) 
with n = 1 for the level 2 data x\, and p\. We continue 
in this fashion to obtain the data at levels 3, 4, etc. 

Strictly speaking, neither of the options outlined above 
is an initial value problem. With the first option, the 
freely specifiable data p\ are split between the initial 
time node and the first time zone. With the second op- 
tion, the data x° a and x\ are split between time nodes 
and 1. Apart from this slight misuse of the word "ini- 
tial" , we see that it is fairly trivial to reinterpret the vari- 
ational integrator Eqs. I110[> as an initial value problem. 
With higher order discretizations, this reinterpretation is 
not so simple [l3|. 



IV. SYMPLECTIC FORM, NOETHER'S 
THEOREM AND THE MIDPOINT RULE 

In this section we show that the variational integrator 
(fTUf) is symplectic, and that Noether's theorem applies. 
These results are derived in mathematically precise lan- 
guage for the Lagrangian formulation of mechanics by 
Marsden et al. [HJ, Hj- In the process of developing 
these results, we show that the VI equations can be ex- 
pressed in terms of the node-centered momentum. The 
discrete equations are equivalent to the midpoint rule ap- 
plied to Hamilton's equations. 

Consider first the continuous Hamiltonian system de- 
fined by the action QJ. The canonical two-form is de- 
fined by 



w = dp a A dx a 



(11) 
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where d is the exterior derivative and A is the exterior 
product. Hamiltonian systems are symplectic, mean- 
ing that the form to is invariant under time evolution. 
We can derive this result by noting that, for a solution 
of the classical equations of motion the variation 
of the action reduces to the endpoint terms in Eq J5J). 
Let S{x" , t"; x' , t') denote the action evaluated along the 
classical history with endpoint data x a (t') = x' a and 
x a (t") = x" a . We see that dS(x", t"; x', t')/dx% = p a (t"), 
and dS(x",t";x',t')/dx' a = -p a (t')- Then the exterior 
derivative of the action is given by 



dS(x",t";x',t') = Pa dx a 



t" 



(12) 



where p a is the canonical momentum evaluated along the 
classical path. The identity dd S(x" ,t" ; x' ,t') — shows 
t" 



that dp a A dx a 



vanishes, so that to is constant in time. 



Now turn to the discrete system defined by the action 
(JSJ. Let us define the coefficient of 5x^ that appears in 
6S [Eq. ©] as <->"Pf , where 



' n 



At f dH 

T \dx~ a 



(13) 



Similarly, we can define the coefficient of 8x° a as <+) "P"j 
where 



(+) 



■p 

• n 



„+i At ( dH 
^ + 2 \dx a 



n+l 



(14) 



The VI equations (|10fl define an evolution in phase space. 
Obviously this evolution can be extended to values of 
n beyond nodes and N. Likewise, we can apply our 
definitions of ( ~ )- P" and <+), P™ for all integer n. Now 
observe that the VI equations imply that w V n a — ( ~ )- P" 
vanishes when the extended equations of motion hold. 
Thus, we can drop the superscripts (+) and (— ) and 
denote both ^V n a and ^V n a by V%. 

Let S(x , t ; x°, t°) denote the value of the discrete 
action (J5J for a solution of the VI equations of motion 
with endpoint data x^ at t N and x° a at t°. The variation 
in Eq. JSJ shows that, when the (extended) equations of 
motion hold. 



ds{x N ,t N -x^f) = vy x n a 



TV 



n=0 



(15) 



This is the analog of Eq. I|12l) above. Taking the exterior 
derivative of this expression we find 



= dV"Adxi 



N 



(16) 



Thus, the discrete action naturally defines a symplectic 
two-form 



uj = dV™Adx" 



(17) 



that is conserved under the phase space evolution defined 
by the VI equations of motion. 



In the analysis above we defined 



P- 



n+l 
PS. 



2 \dx a 



At 



/ dH 
2 \dx^ 



n+l 



(18) 



The two expressions for T 3 ™ are equivalent when the equa- 
tions of motion hold. A short calculation using Eq. (|10b(l 
shows that 



V. 



n+l 



n+l 



(19) 



Therefore we see that, when the equations of motion hold, 
■p™ can be identified as the node momentum p™. 

The equation of motion for T>™ can be derived by com- 
puting AT™ an d using the VI equation ijlObfl . Along with 
Eq. l|lUa|) . we have 



Ax r > +1 (dH 

At ~ \dp~ a 

AV2+ 1 _ ( OH 

At ~ ~ \dx~ a 



n+l 



n+l 



(20a) 
(20b) 



This is perhaps the most elegant form of the VI equa- 
tions. They are simply Hamilton's equations discretized 
with the midpoint rule. Their interpretation as an initial 
value problem is straightforward: given data x a and V® 
at the initial time, we solve the equations with n = for 

x\, V\. Repeat to find data at nodes n = 2, 3, Recall 

that the derivatives of H that appear on the right-hand 
sides of Eqs. I|20() are zone centered functions. Thus, they 
are evaluated at Vtt— xir—, and £2±i. 

Noether's theorem states that symmetries give rise to 
conserved "charges". We now show that when the dis- 
crete action is invariant under a symmetry transforma- 
tion, there exists a charge that is exactly conserved by 
the VI equations. 

Consider first the continuum case. Let x a — > X%(x) be 
a one-parameter family of transformations that leave the 
action, expressed in Lagrangian form, unchanged. Since 
we are working with the Hamiltonian formalism, let us 
extend this family of configuration space transformations 
to a family of point canonical transformations: 



Pa 



a ~ Pb dx a 



(21a) 
(21b) 



Here, it is assumed that a — coincides with the 
identity transformation. By differentiating the relation 
x a = X- a (X a (x)) we see that P£ = p a x a so the 
transformation (|21|l is indeed canonical. 

By assumption the action Q is unchanged by the 
transformations (|21|l . so we have S[p, x] — SIP' 7 , X a ] for 
all a. It follows that the derivative of S[P a , X a ] with re- 
spect to a must vanish. On the other hand, the endpoint 
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terms in the general variation of the action J5J imply that, 
if X%, P° satisfy the equations of motion when a = 0, 
then 



dS[P°,X e 
da 



dXi 



Pa - 



da 



t" 



(22) 



at a = 0. Because the left-hand side of this equation 
vanishes, we see that the charge 



Q = Pa 



dX c n 



da 



(23) 



is conserved in time for the classical motion of the system. 

Now turn to the discrete case. Let us assume that the 
discrete action (J5J is unchanged when the variables a;™, p™ 
are transformed by Eqs (|21|l for each value of n. Then, as 
in the continuum case, the derivative of S[P a ,X a ] with 
respect to a vanishes. The general variation of the action 
© implies that, if (X°) n , (P£) n satisfy the equations of 
motion at a = 0. then 



dS[P a 1 X a ] 
da 



= VI 



Ax a aY 



da 



N 



n=0 



at a = 0. Here, I have used the definitions (|18|l for P™. 
Since the left-hand side of this relationship vanishes, we 
find that the charge 



Q = Vl 



da 



(25) 



is conserved (independent of n) by the VI Eqs. I|10|) or 



V. ENERGY CONSERVATION 

One of the key characteristics of variational integra- 
tors that makes them interesting and important is their 
behavior with respect to energy. If the Hamiltonian has 
no explicit time t dependence then the energy is con- 
served in the continuum theory; see Eq. QJ. Variational 
integrators do not conserve energy exactly, but typically 
the energy error does not grow as the evolution time in- 
creases. To be precise, in this section I show that the VI 
equations (|2*U|) exactly conserve the value of a phase space 
function TL that differs from the Hamiltonian H by terms 
of order 0(At 2 ). The coefficient of the 0(At 2 ) difference 
is a phase space function that remains bounded at least 
as long as the solution trajectory is bounded in phase 
space. It follows that the value of energy predicted by 
the VI equations will be "close" to the exact value, where 
"close" means that the error is of order 0(At 2 ) with a 
coefficient that does not exhibit unbounded growth in 
time. 

Let us begin the analysis by considering the continuum 
evolution for a system with time-independent Hamilto- 
nian H.(p,x). This system is described by the action 



S[p, x] = £ dt [p a ±a - H(p, x)\ 



(26) 



with variation 



5S = eom's + Pa&Xa 



(27) 



The terms listed as "eom's" are the terms that yield 
Hamilton's equations of motion. Let S(x" ,t";x',t') de- 
note the action l|2*6"|l evaluated at the solution of Hamil- 
ton's equations with endpoint data x' a at t' and x" a at 
t". The variation Eq. J^J) implies that S(x",t";x',t') 
satisfies 



dS(x'\t";x',t') 
dx'l 

dS(x" 1 t , ';x',t') 
dx< 



-P' , 



(28a) 
(28b) 



where p' a = p a (t') andp" = p a (t"). These equations show 
that —S(x" ) t"; x' , t') is a Type 1 generating function for 
a canonical transformation from "old" coordinates and 
momenta x' a , p' a to "new" coordinates and momenta x", 
(24) Pa -24|. We also know that, starting from the initial data 
x' a , p' a , the classical trajectory generated by the Hamilto- 
nian Tt(p, x) passes through the phase space point x", p"- 
Thus, —S(x",t";x',t') is a Type 1 generating function 
that generates a canonical transformation representing 
the time evolution of the system from t' to t" . 

Generating functions of different type are related by 
functions of the old and new coordinates and momenta. 
We can define a new generating function H by 



H 



_p'L+p' a <-< S(x",t";x>,t>) 



t" - t' 



t" - t' 



(29) 



Equations l|28ll can be written as dS(x" ,t";x',t') = 
p'adx'a — p' a dx' a . From this result it is straightforward 



to show that the exterior derivative of H is given by 



dH 



A.Xg 

At 



dp a 



A Pa 
At 



dx„ 



(30) 



-Pa, 



where p a = {p^+p' a )/2, x a = « + <)/2, Ap a =p" a - 
Ax a = x" a — x' a , and At = t" — t' . Thus, H can be viewed 
as a function of p a and x a . The canonical transformation 
that represents the classical evolution from t' to t" is 
written in terms of the new generating function H = 
H(p, x) as 



dH(p,x) Ax a 

dpa 
dH(p,x) 



dx a 



At ' 

Ap a 

At 



(31a) 
(31b) 



These are precisely the VI equations (|2Up. the midpoint 
rule, with some simple changes of notation. 

The analysis above shows that the VI equations 12U|) 
can be viewed as the generating function equations for a 
canonical transformation from old coordinates and mo- 
menta x™, T 5 ™ to new coordinates and momenta x™ + , 



The generating function is H(V^1, 



„n+l 



At); it 
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is helpful at this point in the analysis to consider H as 
dependent on the timestep At as well as the coordinates 
and momenta. The canonical transformation generated 
by H defines a mapping of phase space that coincides 
with the exact time evolution described by the Hamil- 
tonian Ti(V,x). The relationship between H and Ti is 
given by 



11+ 



iAx" +1 



At 

S(x n+1 ,t n+1 ;x n ,t n ) 
At 



(32) 



This is Eq. (|29[1 with appropriate changes in notation. 
The function S(x n+1 , t n+1 ; x n , t") is the continuum ac- 
tion l|2t)|) evaluated at the solution of the equations of 
motion with endpoint data at t n and at t n+1 . 

Analogous to Eqs. I|28[) . we have the relations 



•pii+i 

' a 

' a 



dS(x 



n+l £ti+l. x n 



dx2 +1 

dS(x n+1 ,t n+1 ;x n ,t n ) 



da 



(33a) 
(33b) 



that define the momenta at the endpoints. 

The discrete evolution defined by the VI equations with 
Hamiltonian H coincides with the exact continuum evo- 
lution defined by Hamilton's equations with Hamiltonian 
Ti. Since the exact evolution conserves Ti, it follows that 
the VI equations conserve Ti. We now show by explicit 
calculation that H and Ti differ by terms of order 0(At 2 ). 

In order to evaluate the action S along the classical 
solution between t n and t n+1 , we first expand the solution 
x a {t), Pa(t) in a series in t with coefficients that depend 
on X7i — and — = Ve — • The calculation is simplified 
by writing Hamilton's equations Q as 



£a — U ab Ti b 



(34) 



where £ a denotes the set of canonical variables {x a ,p a } 
and LOab is the matrix 



-I 

1 



(35) 



Here and below, subscripts on Tt denote derivatives; for 
example, Ti b = dTi/d^ b . The solution is 

+Uaa> [H a < +H albc LJwH b ,LU ccl H d At 2 /8] (t - <2±i) 

+\u aa .H a , b ujwH b . [(* - t^-f - {At n+1 ) 2 /A\ 



24 



■Jaa' [Tta'bcWbb'Ttb'Wcc'Ttc' + Tta'b^Jbb'Ttb'c^cc'TLc'] 



. [4( t _ tli±l) 3 - 3(< - t^)(At n+1 ) 2 
-0(At 4 ) , 



(36) 



where all derivatives of Ti are evaluated at £a 



n+l 



The Type I generating function -S(x n+1 ,t n+1 ;x n ,t n ) 

n+l 

is written as a function of £s — by inserting the solution 
(|36|l into the action (|26|l . with initial and final times t n 
and t n+1 . The new generating function H is then found 
from Eq. 132|) . with the result 



H = TL 



I 

24 



TLabtOaa'TLa'OJbw'HvAt 2 + 0(At 4 



(37) 



Clearly, this formal expansion for H in terms of Tt can 
be inverted to yield 



I 2 
TL — H- — H ab uj aal H a >Lu bb > H b > At 



0(At 4 



(38) 



This is the desired relationship between the phase space 
functions Ti and H. 

With the solution £ a (t) expanded to terms of order At 3 
in Eq. Q36[) . the evaluation of Eq. l-'il'l) yields H through 
terms of order At 2 . However, a simple argument can 
be given to show that the terms of order At 3 , and in 
fact all terms proportional to odd powers of At, must 
vanish. Consider Eq. (|3*2")l . but let the data t n , x™ and 
t n+1 , exchange roles. The data must be exchanged 
in the definitions i.VMi as well; this yields 



->n+l _ 



dS{x n ,t n ;x n+1 ,t n+1 ) 
dxl 

dS(x n ,t n ;x n+1 ,t n+1 ) 



dxl 



n+l 



(39a) 
(39b) 



Now, the function S(x n ,t n ; x n+1 , t n+1 ) is just the action 
evaluated at the solution of Hamilton's equations with 
endpoint data x a {t n ) = x n a and x a (t n+1 ) = x™ +1 . It 
differs from S(x n+1 , t n+1 ; x n , t n ) only because the limits 
of integration are reversed. Hence, we have 

S{x n ,t n ;x n+1 ,t n+l ) = -S(x n+1 ,t n+1 ;x n ,t n ) , (40) 

and we find that the definitions H39fl & re identical to 
Eqs. ■ It follows that the right-hand side of Eq. J32J| 
is unchanged when we exchange the endpoint data. 
Equating the left-hand sides leads to 

H(V^, a£±i, At) = H CP2±i x^+i, -At) . (41) 

Therefore, H is an even function of At, and it's expansion 
(|3*7ll in terms of Ti does not contain odd powers of At. 

To summarize, the VI equations exactly conserve Ti 
and the Hamiltonian H differs from Ti by terms of or- 
der At 2 . The coefficient of the 0(At 2 ) and higher or- 
der terms are constructed from derivatives of H. As 
long as the motion in phase space remains bounded, and 
the Hamiltonian and its derivatives are nonsingular, then 
these coefficients will remain bounded. It follows that H 
will remain "close" to Ti, which is constant, for all time. 
If the motion in phase space does not remain bounded, 
it does not necessarily follow that the coefficient of the 
0(At 2 ) will grow in time. In this situation the results 
depend on the details of the Hamiltonian. 



8 



The energy behavior of the VI equations is quite differ- 
ent from the behavior exhibited by many numerical inte- 
grators. For example, second-order Runge-Kutta (RK2) 
typically exhibits errors in H of order At 2 on short time 
scales and a drift in the value of H of order At 3 on long 
time scales. Fourth order Runge-Kutta (RK4) exhibits 
errors in H of order At 4 on short time scales and a drift 
in H of order At 5 on long time scales. For both RK2 
and RK4, the energy error becomes unboundedly large 
as time increases. We will see examples of these behav- 
iors in the next section. 



VI. EXAMPLES 

The examples in this section show the results obtained 
from numerical integration of simple Hamiltonian sys- 
tems using the VI equations l|20[l and standard second 
and fourth order Runge-Kutta. Issues of efficiency are 
ignored. Clearly the midpoint rule, being implicit, is nu- 
merically more expensive to solve than explicit integra- 
tion schemes. However, the aim of this paper is to investi- 
gate the properties of variational-symplectic integrators 
without concern for details of implementation. This is be- 
cause, ultimately, we would like to apply these methods 
to theories like general relativity for which standard inte- 
gration techniques are inadequate. The main goal, then, 
is to find a numerical algorithm that works — efficiency is 
at most a secondary concern. 

One can solve the implicit VI equations using a 
Newton-Raphson method. But in practice it is much 
simpler and more reliable to iterate the equations until 
the answer is unchanged to a prescribed level of accu- 
racy. Thus, given x™ and P™, we begin the first iteration 
with the approximation x" +1 w x™, V™ +1 ~ "Pa- This is 
inserted into the right-hand sides of the VI equations to 
yield improved approximations for x" +1 and ■ The 
whole process is repeated until the desired level of accu- 
racy is achieved. For the higher resolution runs presented 
below, about 5 iterations were needed to reach a solution 
that was accurate to 1 part in 10 13 . For the lower res- 
olution runs, around 15 iterations were needed to reach 
this same level of accuracy. 



A. Coupled harmonic oscillators 

Our first example is the system of coupled harmonic 
oscillators with Hamiltonian 

H = \ (pI + pI) + \{x 2 + y 2 ) + \{* - vf ■ (42) 

The graph in Fig. [21 shows the amplitude of one of the 
oscillators, function of time. The behavior ex- 

hibited is rather complicated. The two solid curves show 
the results of numerical integration with the VI equa- 
tions (|20[1 and second order Runge-Kutta (RK2), both 




FIG. 2: The amplitude x for the coupled harmonic oscillator 
as a function of time. The VI simulation produces the solid 
curve that tracks the "exact" solution (dashed curve) fairly 
closely. The other solid curve is obtained from RK2. 

using a timestep of At = 0.1. The dashed curve is ob- 
tained from a fourth order Runge-Kutta integrator with 
timestep At — 0.01. Over the short time scale (t < 50) 
shown in the figure, the dashed curve can be taken as the 
"exact" solution. Compared to RK2, VI does a visibly 
better job of tracking the solution. 

The initial data chosen for the coupled oscillator is 
x = 1, y = p x — p y — 0, so the exact solution has 
energy H = 0.7. Figure |3] shows the error in energy 
for VI at two resolutions. The solid curve shows the 
error for At = 0.01 while the dashed curve shows the 
error divided by 100 for At = 0.1. The close agreement 
between the amplitudes of these two curves shows that 
the energy error is second order in the timestep. The key 
observation is that the energy error does not grow in time, 
even for the simulation with a relatively low resolution 
of At = 0.1. 

The solid curve in Fig. 0] shows the energy error for 
RK2 at a resolution of At — 0.01. The dashed curve 
is the energy error for RK2 with resolution At = 0.02, 
divided by 8. Note that the two curves in this figure 
coincide on long time scales (t > 5). This shows that 
the drift in energy is order At 3 . The short time scale 
errors are 0(At 2 ), so the "wiggles" in the low resolution 
simulation (having been divided by 8) are approximately 
half the size of the wiggles seen in the high resolution run. 
For this particular system, and this particular choice of 
initial data, the growth rate of the energy error with RK2 
is about 2.5At 3 energy units per time unit. 

Qualitatively similar results are found for RK4. In 
Fig.[5]the solid and dashed curves are obtained from sim- 
ulations with timesteps At = 0.01 and 0.02, respectively. 
The errors for the low resolution case have been divided 
by 32. We see that the long time scale drift in energy is 
C(At 5 ), while the short time scale "wiggles" are 0(At 4 ). 
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FIG. 3: Energy error for the coupled harmonic oscillator for 
VI. The solid curve has timestep At — 0.01. The dashed curve 
shows the error divided by 100 for At = 0.1. The results are 
displayed for the first and last ~ 12 time units; the total run 
time was i = lx 10 6 . 
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FIG. 4: Energy error for the coupled harmonic oscillator for 
RK2. The solid curve has timestep At — 0.01. The dashed 
curve shows the error divided by 8 for At = 0.02. 



For this simulation the growth rate of the energy error is 
about —1.1 At 5 energy units per time unit. 

The value of energy H obtained from VI is nearly 
constant because the VI equations exactly conserve the 
nearby Hamiltonian TL. This can be confirmed by com- 
puting the first two terms in the expansion for TL given 
in Eq. (|38Jl . For the coupled harmonic oscillator with 
timestep At = 0.01, the two-term approximation for TL 
remains nearly constant with variations at the level of 
10~ 9 . With timestep At = 0.1, the approximation for 
TL remains nearly constant with variations at the level of 
10~ 5 . These variations are just what we expect given the 



FIG. 5: Energy error for the coupled harmonic oscillator for 
RK4. The solid curve has timestep At = 0.01. The dashed 
curve shows the error divided by 32 for At = 0.02. 



fact that, according to Eq. |JSSJ, the terms omitted in the 
approximation for TL are order 0(At 4 ). 



B. Simple pendulum 

For our next example, consider the simple pendulum 
with Hamiltonian 



H 



1 



p 2 — cos(x) 



(43) 



where x denotes the angle from the vertical and p is the 
angular momentum. Figure [5] shows a portion of the 
phase space for the system. We consider a family of initial 
data points clustered about x = tt/2, p = 0. Specifically, 
the initial data are given by 

x = tt/2 + 0.002 cos(0) , p = 0.002 sin(0) (44) 

for < 9 < 2tt. These points form a "circle" in phase 
space. The boxes in Fig.Elmark the points 9 — 0, n/2, n, 
and 37r/2. The initial data are evolved with VI and with 
RK2, both at low resolution (timestep At = 0.1) and 
high resolution (timestep At = 0.05). The run time is 
74.1 time units, which is just under 10 oscillation periods. 
The initial data cycles around the phase space diagram 
in a clockwise direction. Figure shows the end result 
of this evolution for the two integrators at low and high 
resolutions, as well as the "exact" solution obtained from 
RK4 with a very small timestep. The dashed curves show 
the constant energy contours with energies determined by 
the initial data shown as boxes. 

Qualitatively, we see that both VI and RK2 schemes 
are second order accurate. That is, the errors in x and 
p are reduced by a factor of about 4 when the resolution 
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angle x 

FIG. 6: Phase space diagram for the pendulum. The initial data occupy a "circle" around x — it/2, p = 0, and are evolved for 
just under 10 oscillation periods. Final data is shown for VI and RK2, for low and high resolutions. 



is doubled. But the character of that error is very differ- 
ent. The VI evolution stays close to the constant energy 
contours, and the phase space errors lie almost entirely 
in the H = constant subspace. The RK2 integrator does 
not respect conservation of energy, and over time the sys- 
tem point in phase space spirals outward with increasing 
energy. After about 9000 time units the simulation with 
RK2 and At = 0.1 predicts that the pendulum will gain 
enough energy to circle around completely, rather than 
oscillate. 

Recall that the midpoint rule is a symplectic integra- 
tor, that is, the symplectic form l|17|) is preserved in time. 
It follows that the volume of phase space bounded by the 
initial data "circle" in Fig. H3 is constant under the dis- 
crete evolution defined by the variational integrator. The 
standard second order Runge-Kutta scheme is not sym- 
plectic, and does not preserve phase space volume. In 
Fig.Elit is not possible to tell, simply by looking, whether 
or not the initial phase space volume is conserved by the 
VI scheme, or changed by the RK2 scheme. A more 
involved numerical test would be needed to verify the 
expected results. 



C. Unbounded motion in one dimension 

The VI equations conserve the phase space function TL 
exactly, but the energy H might not remain close to TL 
if the motion of the system is unbounded. Consider the 
Hamiltonian for a particle moving in a one dimensional 
potential, H — p 2 /2 + V(x). In this case Eq. I|38|l gives 

At 2 

H-TL= — [p 2 V" + (V) 2 ] + G(At 4 ) , (45) 



where prime denotes d/dx. The time derivative of this 
difference is d(H - TL)/dt = p 3 V"'At 2 /24: plus terms of 
higher order in At. We see that H — TL, and therefore 
also H, will grow in time if p 3 V" remains finite and does 
not change sign. 

A nice example of this unbounded behavior is ob- 
tained with the potential V(x) = — x 6 / 5 . In this case 
the particle motion at late times is given approximately 
by i ~ (2V2t/5) 5/2 , p ~ V2(2v^t/5) 3 / 2 . Equation 
(I45|l shows that the energy grows linearly with time, 
H ~H+(2V2At 2 /l25)t. Figure □ confirms that for this 




5 10 15 20 

time 



FIG. 7: Energy as a function of time for a particle in a one- 
dimensional potential, V{x) = —x 6 ^ s , obtained with VI. As 
expected, the error in energy grows linearly with time. 

system, the variational integrator exhibits linear growth 
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in the energy. The initial data used in this simulation 
was x — 0.1, p — 0, with timestep At = 0.1. The energy 
error obtained from RK2 is almost identical to the result 
shown in Fig. [7| for VI. 



D. Orbital motion 

Our final example is motion in a gravitational (or elec- 
tric) field described by a central 1/r potential. The 
Hamiltonian is defined by 



H 



1 



i 



(46) 



This system is symmetric under rotations in the x—y 
plane. The conserved Noether charge associated with ro- 
tational symmetry is angular momentum, J — xpy — yPx- 
The initial data for this simulation is x — 1.0, p x = 0.0, 
y = 0.0, p y = 1.2. The resulting orbital motion is an 
ellipse with eccentricity ^0.5 and period ~ 15. 

Figure[S]shows the angular momentum as a function of 
time for RK2 and VI with timestep At = 0.25. With the 
variational integrator the angular momentum is exactly 
conserved (to machine accuracy) and J retains its initial 
value of 1.2 throughout the simulation. With RK2, the 
angular momentum exhibits short timescale fluctuations 
and a longer timescale drift. A more complete analy- 
sis shows that the short timescale errors are order At 2 , 
whereas the drift in J is order At 3 . Qualitatively simi- 
lar results are obtained for RK4. In that case, the short 



timescale errors in J are order At 4 , and the drift is order 
At 5 . 
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FIG. 8: Angular momentum as a function of time for motion 
in a central potential. The solid curve is obtained from RK2. 
The constant, dashed line is obtained with the variational 
integrator. 
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